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Abstract 

At the preliminary design stage of a wing structure, an efficient simulation, one needing little 
computation but yielding adequately accurate results for various response quantities, is essential in 
the search of optimal design in a vast design space. In the present paper, methods of using 
sensitivities up to 2 nd order, and direct application of neural networks are explored. The example 
problem is how to decide the natural frequencies of a wing given the shape variables of the 
structure. It is shown that when sensitivities cannot be obtained analytically, the finite difference 
approach is usually more reliable than a semi-analytical approach provided an appropriate step size 
is used. The use of second order sensitivities is proved of being able to yield much better results 
than the case where only the first order sensitivities are used. When neural networks are trained to 
relate the wing natural frequencies to the shape variables, a negligible computation effort is needed 
to accurately determine the natural frequencies of a new design. 
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Introduction 


The modal response of wing structures is very important for assessing their dynamic response 
including dynamic aeroelastic instabilities. Moreover, in a recent study ' an efficient structural 
optimization approach was developed using structural modes to represent the static aeroelastic wing 
response (both displacements and stresses). 

Sensitivity techniques are frequently used in structural design practices for searching the 
optimal solutions near a baseline design 2 3 . The design parameters for wing structure include 
sizing-type variables (skin thickness, spar or rib sectional area etc.), shape variables (the plan 
surface dimensions and ratios), and topological variables (total spar or rib number, wing topology 
arrangements etc.). Sensitivities to the shape variables are extremely important because of the 
nonlinear dependence of stiffness and mass terms on the shape design variables as compared to the 
linear dependence on the sizing-type design variables. 

Kapania and coworkers have addressed the first order shape sensitivities of the modal response, 
divergence and flutter speed, and divergence dynamic pressure of laminated, box-wing or general 
trapezoidal wing composed of skins, spars and ribs using various approaches of determining the 
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response sensitivities 

In the present paper, the natural frequencies of general trapezoidal wing structures are to be 
approximated using shape sensitivities up to the 2 nd order, and different approaches of computing 
the derivatives are investigated. The baseline design and shape sensitivities are calculated based on 
an equivalent plate-model analysis (EPA) method developed by Kapania and Liu 11 . For 
comparison, an efficient method that employs the artificial neural networks to relate the natural 
frequencies of a wing to its shape variables is also established. An example of a 4 3 full factorial 
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experimental design, i.e., 4 levels in 3 variables, is treated by these methods to display their 
respective merits. 

Shape Sensitivities 


For a trapezoidal wing, there are four major independent shape variables: 1 ) the sweep angle A , 
2) the aspect ratio a , 3) the taper ratio I , and 4) the plan area A . All the other dimensions of the 
wing plate configuration can be calculated using these parameters as follows: 

s = y[aA,a = 2i$/a(\ + T),b = 2s/a(\+T) ( 1 ) 

where s is the length of semi-span, a and b are the chord-length at wing tip and root 
respectively, as shown in Fig. 1 . 

The sensitivities for the design parameters at a baseline design point indicate trends in the 
response of the baseline design if the parameters are perturbed. Usually, only the first order 
derivatives are used: 


f(x',x\---,x n ) = + -x' 0 ) 

i = 1 &X 


(2) 


where 


a/o _ jjf 

dx l dx‘ 


is the sensitivity at the baseline point with respect to the / -th design 
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parameter. For a more accurate approximation, we can use higher-order derivatives in the Taylor’s 
expression: 
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where besides the first order derivatives, second order derivatives . (i, j = are also 


used. 

Equivalent Plate Analysis (EPA) of Trapezoidal Wing Structures 

In Kapania and Liu 1 1 , a general method is developed to analyze trapezoidal wing structures 
composed of skins, spars and ribs. The method is based on the Reissner-Mindlin model, a First- 
Order Shear Deformation Theory (FSDT). An equivalent plate model is based on the hypothesis 
that the original complex built-up structures behave like a plate, a simplification to reduce the 
computational effort to obtain the free vibration and static responses. Compared with the methods 
presented earlier in Kapania and Lovejoy 12 and Cortial 13 , the formulation in Kapania and Liu 1 1 
entails no limitation with respect to the wing thickness distribution. As shown in Ref. 1 1, the 
method shows a good performance for both static and vibration problems in comparison with the 
FEA. Free vibration and static response are obtained using the Ritz method. The advantages of 
using Ritz method are: (i) the ease of carrying out formulation and computer realization, and (ii) the 
global nature of the solution, important for accurate determination of various stress components. 

Due to its efficiency in determining the natural frequencies and mode shapes of wings, the 
Equivalent Plate Analysis (EPA) mentioned above can be used to investigate the variation of modal 
response, that is, to evaluate the sensitivities of the natural frequencies with respect to trapezoidal 
wing structures shape changes. For determining the response of the baseline design, the EPA can 
be used, or the FEA employing a commercial package such as MSC/NASTRAN can be used for 
better accuracy. 
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A key problem that needs to be addressed before this evaluation can be made is mode tracking. 
The natural frequencies given by an ordinary eigenvalue solver are usually ranked by magnitude 
but not by the modal content. As design variables are perturbed, frequencies drift and mode 
crossing may occur. An algorithm for mode tracking is needed to maintain the correspondence 
between eigenpairs of the baseline and the perturbed design. Several methods for such purpose 
have been given by Eldred et al for self-adjoint 14 and nonself-adjoint 15 eigenvalue problems. 

In the present study, a simple yet effective method is used. In this method, any ordinary 
eigenvalue solver can be used, and the modes of the baseline structure are chosen as the 
benchmarks. By using the modal assurance criterion (MAC) defined as 


MAC,, 




( 4 ) 


where \<j ) j } and {0, } are the eigenvector for the perturbed and the baseline design respectively, if 


MAC Jt = max (MAC ;i ) , we say that the j -th mode of the perturbed design corresponds to the i -th 
mode of the baseline structure. 


Approaches to Sensitivity Evaluation 

There can be three kinds of approach for obtaining sensitivity derivatives: the finite difference 
approach, the analytical approach, and the semi-analytical approaches. The finite difference 
approach is very simple to formulate and implement, but is numerically inefficient and is sensitive 
to the step-size used. A too-large step size usually causes significant truncation errors and a too- 
small step size may lead to large round-off errors. As a result, the more elegant and accurate 
analytical approach is used if it does not involve complex mathematical derivation. But for most 
practical problems, the derivation of analytical derivatives is too formidable to handle manually. 
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The basic idea behind the Automatic Differentiation (AD) is to let a computer to perform such 
extensive tasks. The advantage of AD is to avoid truncation errors. The method has found 
applications in sensitivity evaluation 9,10 . For the basic theory of AD one can consult Ref. 16, and 
for the state-of-the-art of AD one can refer to Ref. 17. If an approach uses both analytical and 
finite-difference solutions to obtain the derivative, then it can be called a semi-analytical one. 

The finite difference approaches can be constructed using the following formulas: 


/ , w = / ( x + Ax ) -/(x-Ar) +0 ( Al ,) 

2Ax 


/"(x) = f(x + Ai> - 2/( f> + f{x - W +Q (\ ,■) 

Ax- 


fS) 

( 6 ) 


where 


Ax = £x (7) 

in which £ is the relative step size, but herein it is simply called the step size. Eq. (5) can be 

d 2 f 

applied twice for evaluating the mixed second order derivatives such as — — 1 — (i * j ) . 

3x'3x J 

The analytical approaches for shape sensitivities of modal response can be based on the 
following equations 
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here 6 is the shape variable, A. and {0, } are the i -th eigenvalue and eigenvetor, and [<j ) t } is mass- 
normalized such that {0, } r [M]{0, } = I . Equations (8) and (9) were first derived by Wittrick 18 and 
Fox and Kapoor 19 respectively. One can find more on this topic in Ref. 20. 


The major difficulty of applying Eqs. (8) and (9) lies in the calculation of 


am 

36 


and 


3 [M ] 
36 


For instance, consider 


am 

a6 


. According to Ref. 1 1, the stiffness matrix [A'] is formulated as an 


integral 


[K] = JJJ [C] T [7] r [D][T][C]dV = J' J‘ [C] T ( f* [T] t [D][T]\ j\dz \c]d^dr] 

v ~ \ z ' J (11) 

= j\j\[C] T [G][C]d&7] 


where only the inner part [G] = J 2 [7] r [D][r]|y|dz is a function of the shape variables, and the 
Gaussian quadrature is used to obtain the integration on £ and 77 . Therefore, 


in which 


3 e { do 7 

a[G] 


( 12 ) 


36 


can either be determined analytically or numerically. 


Often people make use of the advantages of both the finite difference and analytical approaches 
in different stages of obtaining some complicated sensitivities. While trying to use the analytical 
approach as much as possible, in other parts of the process the finite difference is used, as in the 
cases of Refs. 8 and 10. This kind of approach is usually called semi-analytical. 

In summary, there are three approaches to calculate the first order modal sensitivities: 
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(i) analytical approach : Eqs.(8)~(10) are used, and 1 and — 4~~ ~ are determined 


36 


36 


analytically. 


(ii) semi-analytical approach : Sensitivities and j n Eqs.(8)~(10) are determined 


36 


36 


k i i 

numerically, that is, for the case of — — , Eq. (12) is used where— — -is calculated using a finite 

36 36 

difference scheme. 

(iii) finite difference approach: and — ~ - are determined using Eq. (5) directly. 

3 6 36 

For the second order sensitivities, there can still be three approaches as specified above. While 
the formulation for the analytical approaches is becoming more complicated, a scheme as simple as 
Eq. (6) can be used for the finite difference approach. 


Application of Sensitivity Technique (ST) in Multi-variable Optimization 


In a multi-variable case, the following formulation is used instead of Eq. (3): 


r) R 1 

R i (p) = R(p i ) + (p-p i f— + - 

3 Pi 2 


-|2 


iP -Pi)' 


3 Pi 
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(13) 


where p = (v',v 2 ,---,v' , ) r is an arbitrary point in the design space, p t =(v ( l ,v/,”-,v' 1 )' is the i- 


1 ..2 ..n^T 


th node point in the design space, R t ( p ) is the response at p estimated by using the response and 
its sensitivities at p f , R(p t ) is the response at the i -th node point p t , and 


3 
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Once there are enough estimates for the response at p using Eq. (13), a more accurate 


evaluation of response at p can be determined using the following weighting procedure involving 


the so-called exponentially decaying influence function 21 : 

R(p) = ^Wi(p)Ri(p) 


( 14 ) 


where i ranges through the N w design points which are closest to p , and the weight coefficients 
w l ( p ) are determined such that its sum is unity: 


w,(p) = 


exp(-C|/?-/>,|) 

£exp(-C|p-p,.|) 


( 15 ) 


in which C is an empirical constant, and the distance between p and p s is defined as 


p - p ( =J^(v J - v/) 2 . It can be seen that ^w ( (p) = 1 . 


j = i 


Application of Neural Networks (NN) 


Artificial Neural Networks (ANN), or simply Neural Networks (NN) are computational systems 
inspired by the biological brain in their structure, data processing and restoring method, and 
learning ability. More specifically, a neural network is defined as a massively parallel distributed 
processor that has a natural propensity for storing experiential knowledge and making it available 
for future use by resembling the brain in two aspects: (i) Knowledge is acquired by the network 
through a learning process; (ii) Inter-neuron connection strengths known as synaptic weights (or 
simply weights) are used to store the knowledge 22 . 

The NN has the following properties: (i) Many of its kind are universal approximators, in the 
sense that, given a dimension (number of hidden layers and neurons of each layer) large enough, 
any continuous mapping can be realized; Therefore, (ii) a NN provides a general mechanism for 
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building models from data, or gives a general means to set up input-output mapping (iii) The input 
and output relationship of NN can be highly nonlinear; (iv) A NN is parallel in nature, and it can 
make computation fast when executed in a parallel computer, though NN can be simulated 
in ordinary computers in a sequential manner. 

Major steps of utilizing NN include: (i) specifying the topology or the structural parameters 
(number of layers, number of neurons in each layer, etc.) of the NN, (ii) training of the NN, 
corresponding to the learning process of the brain, (iii) simulation, corresponding to the recalling 
function of the brain. 

In the present work MATLAB Neural Network Toolbox was used. Ref. 23 contains a summary 
of applications of NN in structural engineering and details of how to make use of MATLAB Neural 
Network Toolbox. 

Generally speaking, there can be two directions to use NN for the efficient simulation of the 
performances of wing structures as in the following: 

(a) Direct Application 

In this case, the input layer includes all the design variables of interest (for instance, the four 
shape parameters of the wing plan-form). The output layer gives the desired structural 
responses, such as natural frequencies etc. The EPA is being used as the training data generator, 
though if necessary, results obtained using a FEM can also be used as the training data. 
Preparation of training data is very important, and the training algorithm used also greatly 
impacts the training process 23 . Caution must be exerted in specifying the network parameters 
and training criterion so that the results of the trained network would not oscillate around the 
training data. The direct application is what we do in this paper. 
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(b) Indirect Application 

In this approach, the aim is to develop a way of incorporating NN into the application of the 
equivalent plate model analysis (EPA) of complex wing structures, other than just making use 
of results generated by EPA as the training data base. Note that in the EPA of a complex wing, 
computational efforts lies mainly in performing the various integrals for the inner-structural 
components of the wing, i.e. the spars and the ribs. If an anisotropic material can be found to 
replace this structure such that the new composite wing has similar global properties as the 
original one, then the EPA can be performed more efficiently. Determining the adequate 
material properties of the equivalent anisotropic material is the major obstacle here. The role of 
NN will be in relating the material properties to all kinds of wing design parameters, when there 
exists enough database for training. Use of NN to determine properties of an equivalent 
anisotropic material that will accurately represent the wing is being studied at present in another 
effort. 

Examples and Discussion 

(a) Results on sensitivity evaluation 

Particulars of the baseline wing structure are as follows: the sweep angle A = 30° , the aspect 
ratio a = 3.5 , the taper ratio i = 0.5 , the plan area A = 5832m 2 . The wing sections are generated 
using the Karman-Trefftz transformation (Ref. 24) and has a thickness-chord ratio of 0. 15 at the 
wing root and 0.06 at the tip. The skin thickness t 0 = . 118m . There are 4 spars and 10 ribs 
distributed uniformly under the skins. Particulars of the spars and ribs are the same: the cap height 
/i, = . 197m , the cap width /, = ,373m , and the web thickness r, = .059m . There is only one 
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material used with mass density p = 2.526 x 10^/6 • sec 2 /in*. Young’s modulus 

E = 1 .025 x 1 0 1 lb/ in 2 , and Poisson's ratio v = 0.3 . The wing is clamped at the root. 

An example of using EPA to calculate the natural frequencies with regard to shape variables 
while tracking modes by evaluating MACs is provided in Fig. 2, where the variation of the natural 
frequencies of the first 10 modes w.r.t. the aspect ratio are shown. It can be seen that for most cases 
the intersection of natural frequencies has been treated well, and only in a few cases the frequency 
variations near the intersection point seem to have a minor problem, probably due to some kind of 
interaction between the two modes. If an eigenvalue solver that can work more accurately with 
repeated eigenvalues is made use of, the situation can be improved. 

The effect of step size on the finite difference approach for sensitivities was investigated for all 
the four shape variables. The case with the taper ratio is shown in Fig. 3. From all the cases, it is 
seen that for the best results for both the 1 st and 2 nd order sensitivities, the step size £ defined in 
Eq. (7), should be between 0.005-0.015 and for fairly accurate results £ can be between 
0.0017-0.045. 


To evaluate 


BIG] 

36 


analytically proved to be formidable except only in some simplified cases. 


In order to compare the sensitivities using the analytical, semi-analytical and finite difference 
approach, a special case of the above baseline wing with a constant thickness was considered so 


that the analytical derivation of 


3[G] 

36 


in Eq. (12) is not formidable. When e is specified as 0.005, 


it is found that for the 1 st order sensitivities to the four shape variables ( A, a, z, and A ) the relative 
difference (averaged for the first 10 modes) between the finite difference and analytical approach is 
0.003%, 0.003%, 0.002% and 0.003% respectively. The relative difference between the semi- 
analytical and analytical approach is 0.14%, 0.04%, 0.02% and 0.01% respectively. Therefore in 
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this case the finite difference approach is more accurate than the semi-analytical one, however both 
the approaches yield quite accurate results. 


For the original baseline wing, since the derivation of the analytical derivatives for 


<nc] 

96 


is 


too formidable, only the comparison of the l sl order sensitivities using the finite difference and the 
semi-analytical approach was made. It is found in this case the sensitivities to the aspect ratio a , 
taper ratio i and plan area A using both approaches are quite close, the average difference for the 
first 10 modes being in the range of 0.5~1.4%. As an example, the 2 nd natural frequency w.r.t. A is 
shown in Fig. 4, where it can be seen that the l sl order sensitivities using the finite difference and 
the semi-analytical approach almost coincide with each other. On the other hand, sensitivities to the 
sweep angle A using the two approaches have had some quite large relative differences especially 
for modes whose sensitivity to A is small. One such example, the 3 rd natural frequency w.r.t. A , is 
shown in Fig. 5, where attention should be paid to the scale for the vertical coordinate to see how 
small the sensitivity to A really is. 

It is observed in Fig. 5 that, as in the case of the constant-thickness wing, the finite difference 
approach has a better performance than the semi-analytical one. In fact, in some extreme cases, the 
linear approximation using the first order sensitivity obtained using the semi-analytical approach is 
not at all tangent to the actual variation at the baseline point. This is not the case for that using the 
finite difference approach, if the step size chosen is not too large. Moreover, the computation 
efforts for both the approaches are at the same level since in both cases calculation of the stiffness 
and mass matrices at the baseline design and two perturbed designs should be performed. 

It is obvious from observing Fig. 4 and 5 that the approximation using sensitivities up to the 
second order has much improved the results compared with the case where only the first order 
sensitivity is used. Similarly it has been shown in Haftka and Gurdal 25 that, for the stress-ratio in a 
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three-bar truss, the quadratic approximation is much more accurate than the linear one. Also it can 
be seen that the second order sensitivities using the finite difference scheme of Eq. (6) are fairly 
accurate, at least for the purpose of engineering application. Another advantage of this scheme is 
that it shares the perturbation data with the first order sensitivity scheme Eq. (5), therefore its 
evaluation has no increase in the computational effort at all. 

Using the finite difference approach based on Eq. (5) the mixed second order sensitivities 


d 2 /o , 


dx'dx J 


( i * j ) can be readily determined. As an example, the mixed second order sensitivity on z 


and A for the first five natural frequencies were calculated, and the results are listed as follows: 


0.0099, 0.0153, 0.0353, 0.0494 and 0.0156. 


(b) Application of Sensitivity Technique (ST) and Neural Networks (NN) 

For a trapezoidal wing, there are four major independent shape variables, i.e. the sweep angle 
A , the aspect ratio a , the taper ratio z , and the plan area A . As an example , a 4 3 full factorial 
experimental design with 4 levels in A, a , and z respectively, was used. Particulars of the levels 

of every variable are as follows: A = [0° ,10° ,20° ,30°], a = [1 .0,1 .5,2.0,2.5] , and 


z = [0.3,0.4,0.5,0.6] . The plan area is chosen to be a constant: A = 3500/Vr . The other particulars 
are the same as in (a). 

The natural frequencies of the wing structure at the 64 node points in the design space were 
calculated using EPA, and the 1 st and 2 nd order sensitivities at these points were also determined by 
finite difference using EPA 3 . For each mode, a feed-forward neural network with a structure of 
3x15x10x1, i.e. 3 inputs, 15 neurons in the first hidden layer, 10 neurons in the second hidden 
layer, and 1 output, is trained using the MATLAB NN Toolbox function trainlm that trains feed- 
forward network with the Levenberg-Marquardt algorithm 23 . There are 64 sets of training data, 
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which are non-dimensionalized before the training process. Once the networks are trained, the 
input-output relationships can be readily retrieved by using the function simuff. 

For the application of sensitivity technique, the major task is to evaluate the sensitivities, and to 
generate responses at an arbitrary design point using Eqs. (13) and (14) does not need large amount 
of CPU time. The constant C in Eq. (14) was specified to be 10, and N w = 10 was used. 

Shown in Fig. 6 are the first 6 natural frequencies of 20 randomly chosen wing structures inside 
the box defined in terms of lower and upper bounds on the design variables specified above. From 
the figure it can be seen that both of the results given by NN and ST are in very good agreement 
with the desired values (those given by the EPA) except for a few cases where there are some 
differences. These cases might be caused by the unstable performance of the algorithm used for 
extracting eigenvalues in the EPA near the mode-crossing points, as shall be shown in Fig 7 and 8. 
In order to see the effects of sensitivity order, a randomly chosen path inside the design space box 
is defined as 


v j =v'(\-a i ) + v'a i ,j = 1,2,3' 
v 1 = A,v 2 =a,v 3 = r, 
a i =s n ' ,rij =r y /(l-/v). 


(16) 


where and v{ are lower and upper bounds of variable v j , for instance, v^=0°, vj =30° etc., 


se [0,1] is the range of a shape variable, and r (j = 1,2,3) are randomly determined values between 

0 and 1 . Results of natural frequencies of the first 4 modes for wing structures defined by points 
along a path with n , =0.945, n 2 =8.200, and n 3 =3.203 are shown in Fig. 7, where only the 1 st order 

sensitivities were used, and in Fig. 8, where sensitivities up to the 2 nd order were used. It can be 
seen that when sensitivities up to the 2 nd order are used, results are effectively improved. 
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Generally speaking, neural networks and sensitivity technique can give equally good results, 
and the former uses less time than the latter. But both methods, once the NNs are trained or the 
sensitivities are obtained, are much more efficient than the EPA . For instance, the CPU times 
consumed by the EPA, the sensitivity based method and the NN based method are in the ratio of 
55 : 1 : 0.06 . 

The example used above has only three variables. For design problems with more variables, the 
method of NN and ST can still be applied in general, only at the expense of more computing time. 
We can expect that similar conclusions to those obtained above still apply to these cases. For a 
design problem with very large number of variables, in combination with the NN or ST method, 
methodologies to shrink the design space, such as the reasonable design space approach described 
in Ref. 26, can be used. This can make the search of optimal design easier and at the same time the 
application of NN or ST more accurate, just as the case in Ref. 26 where the response surface 
approximation was used to simulate high-fidelity models. Also for this kind of high dimensionality 
design problems, a full multi-level factorial experimental design is almost impossible to use hence 
the methods of either NN, or ST, or even response surface are hard to apply because the cost would 
be too high. In such a case, an incomplete block statistical experimental design using the D-optimal 
criterion 27,28 can be used, which, with much reduced number of design node points, makes the 
application of NN or ST possible. 

Conclusion 

Modal response of general trapezoidal wing structures was investigated based on an equivalent 
model analysis and sensitivity techniques. The variations of the natural frequencies w.r.t. shape 
design variables need to be coordinated with the baseline mode shapes by mode tracking. The use 


16 



of second order sensitivities proved to be yielding much better results than the case where only first 
order sensitivities are used. Shape sensitivities can be evaluated using analytical, finite difference 
and semi-analytical approaches. The present research shows that when the analytical solution is not 
available, the finite difference approach would be a better choice than the semi-analytical one 
provided the step size is properly specified. 

Neural networks can be trained to relate the natural frequencies of a wing structure to its shape 
variables. In this approach the major efforts are in training the networks. Once the networks are 
trained, there needs an almost negligible computational effort to obtain equally good results for the 
natural frequencies for any given set of the wing shape variables. 
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Fig. 1 Plan configuration of a trapezoidal wing 
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Fig. 4 The 2 nd natural frequency w.r.t. wing plan area 
using 1 st and 2 nd order sensitivities 
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Natural Frequency (rad/sec) 



Fig. 5 The 3 rd natural frequency w.r.t. wing sweep angle 
using 1 st and 2 nd order sensitivities 
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Fig. 6 Comparison of the natural frequencies of the first 6 modes 
for wing structures randomly chosen inside the box 
of design space, as obtained by the NN and ST w.r.t. 
those obtained using a full-fledged EPA 
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Simulated Frequencies (rad/sec) 



Fig. 7 Comparison of the natural frequencies of the first 4 modes 
for wing structures along a path inside 
the box of design space (n , =0.945, n 2 =8.200, n 3 =3.203) 

using only the 1 st order sensitivities 
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Fig. 8 Comparison of the natural frequencies of the first 4 modes 
for wing structures along a path inside 
the box of design space (n , =0.945, n 2 =8.200, n 3 =3.203) 
using sensitivities up to the 2 nd order 
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